This notebook demonstrates the Australian Water Observations from Space (WOFS) algorithm. This water detection algorithm is an improvement over the Landsat QA water flag or the NDWI index for water identification. For more information, visit this website:
Import Dependencies and Connect to the Data Cube ▴¶
In [1]:
importdatacubedc=datacube.Datacube(app='Water_Observations_from_Space')fromdatacube.utilsimportmaskingimportsys,osos.environ['USE_PYGEOS']='0'importdatetimeimportmatplotlib.pyplotaspltimportnumpyasnpimportxarrayasxrimportpandasaspdfromdea_tools.plottingimportrgb,display_mapfromdea_tools.bandindicesimportcalculate_indices### EASI toolssys.path.append(os.path.expanduser('../scripts'))fromceos_utils.data_cube_utilities.clean_maskimportlandsat_clean_mask_invalid,landsat_qa_clean_maskfromeasi_toolsimportEasiDefaultsfromeasi_toolsimportnotebook_utilseasi=EasiDefaults()# Get the default parameters for this system
Successfully found configuration for deployment "eail"
# Select an analysis region (Latitude-Longitude) # Select a time period within the extents of the dataset (Year-Month-Day)# Mombasa, Kenya# latitude = (-4.05, -3.95) # longitude = (39.60, 39.68) # latitude=easi.latitude# longitude=easi.longitudelatitude=(36.28,36.48)longitude=(-114.325,-114.43)# Define Time Range# Landsat-8 time range: 07-Apr-2013 to currenttime_extents=('2021-01-01','2021-12-31')
In [6]:
# The code below renders a map that can be used to view the analysis region.display_map(longitude,latitude)
/env/lib/python3.10/site-packages/dea_tools/plotting.py:313: FutureWarning: This function is deprecated. See: https://pyproj4.github.io/pyproj/stable/gotchas.html#upgrading-to-pyproj-2-from-pyproj-1
all_longitude, all_latitude = transform(Proj(crs), Proj('EPSG:4326'), all_x,
Out[6]:
Make this Notebook Trusted to load map: File -> Trust Notebook
After loading, you will view the Xarray dataset. Notice the dimensions represent the number of pixels in your latitude and longitude dimension as well as the number of time slices (time) in your time series.
Time series output of the Australian Water Observations from Space (WOFS) results. The results show the percent of time that a pixel is classified as water over the entire time series. BLUE = frequent water, RED = infrequent water.
In [13]:
# WOFS is written for Landsat Collection 1, so we need to scale the Collection 2 data to look like collection 1# This is stolen from https://github.com/GeoscienceAustralia/wofs/blob/master/wofs/virtualproduct.pydefscale_usgs_collection2(data):"""These are taken from the Fractional Cover scaling values"""attrs=data.attrsdata=data.apply(scale_and_clip_dataarray,keep_attrs=False,scale_factor=0.275,add_offset=-2000,clip_range=None,valid_range=(0,10000))data.attrs=attrsreturndatadefscale_and_clip_dataarray(dataarray:xr.DataArray,*,scale_factor=1,add_offset=0,clip_range=None,valid_range=None,new_nodata=-999,new_dtype='int16'):orig_attrs=dataarray.attrsnodata=dataarray.attrs['nodata']mask=dataarray.data==nodata# add another mask here for if data > 10000 then also make that nodatadataarray=dataarray*scale_factor+add_offsetifclip_rangeisnotNone:dataarray=dataarray.clip(*clip_range)dataarray=dataarray.astype(new_dtype)dataarray.data[mask]=new_nodataifvalid_rangeisnotNone:valid_min,valid_max=valid_rangedataarray=dataarray.where(dataarray>=valid_min,new_nodata)dataarray=dataarray.where(dataarray<=valid_max,new_nodata)dataarray.attrs=orig_attrsdataarray.attrs['nodata']=new_nodatareturndataarraylandsat_dataset_scaled=scale_usgs_collection2(landsat_dataset)
/env/lib/python3.10/site-packages/rasterio/warp.py:344: NotGeoreferencedWarning: Dataset has no geotransform, gcps, or rpcs. The identity matrix will be returned.
_reproject(
In [15]:
# Apply nan to no_data valuests_water_classification=ts_water_classification.where(ts_water_classification!=-9999).astype(np.float16)# Time series aggregation that ignores nan values. water_classification_percentages=(ts_water_classification.mean(dim=['time'])*100).wofs.rename('water_classification_percentages')
In [16]:
# import color-scheme and set nans (no data) to blackfrommatplotlib.cmimportjet_rjet_r.set_bad('black',1)
# This is where the WOFS time series product is generated. # Areas of RED have experienced little or no water over the time series# Areas of BLUE have experience significant or constant water over the time seriesfigsize=6water_classification_percentages.plot(cmap=jet_r,figsize=(figsize,figsize*img_scale),vmin=0,vmax=100)plt.title("Percent of Samples Classified as Water")plt.axis('off')plt.show()
/env/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: divide by zero encountered in divide
return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: divide by zero encountered in divide
return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: divide by zero encountered in divide
return func(*(_execute_task(a, cache) for a in args))
/env/lib/python3.10/site-packages/dask/core.py:119: RuntimeWarning: divide by zero encountered in divide
return func(*(_execute_task(a, cache) for a in args))
